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ABSTRACT 

We use a multi million particle N-body + SPH simulation to follow the formation of 
^ . a rich galaxy cluster in a ACDM cosmology, with the goal of understanding the origin 

t^J- ' and properties of intracluster stars. The simulation includes gas cooling, star forma- 

Q*^ \ tion, the effects of a uniform ultraviolet background and feedback from supernovae. 

C^> i Halos that host galaxies as faint as = -19.0 are resolved by this simulation, which 

l/"") ' includes 85% of the total galaxy luminosity in a rich cluster. We find that the accumu- 

lation of intracluster light (ICL) is an ongoing process, linked to infall and stripping 
events. The unbound star fraction varies with time between 10% and 22% of the total 
amount of cluster stars, with an overall trend to increase with time. The fraction is 
20% at z = 0, consistent with observations of galaxy clusters. The surface bright- 
O i' ness profile of the cD shows an excess compared to a de Vaucouleur profile near 200 

kpc, which is also consistent with observations. Both massive and small galaxies con- 
tribute to the formation of the ICL, with stars stripped preferentially from the outer, 
lower metallicity, parts of their stellar distributions. Simulated observations of plan- 
, etary nebulae (PNe) show significant substructure in velocity space, tracing separate 

streams of stripped intracluster stars. Despite an unrelaxed distribution, individual 
intracluster PNe might be useful mass tracers if more than 5 fields at a range of radii 
have measured line-of-sight velocities, where an accurate mass calculation depends 
jJJ ■ more on the number of fields than the number of PNe measured per field. However, 

the orbits of IC stars are more anisotropic than those of galaxies or dark matter, 
which leads to a systematic underestimate of cluster mass relative to that calculated 
with galaxies, if not accounted for in dynamical models. Overall, the properties of ICL 
formed in a hierarchical scenario are in good agreement with current observations, 
supporting a model where ICL originates from the dynamical evolution of galaxies in 
dense environments. ICL should thus be ubiquitous in galaxy clusters. 
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1 INTRODUCTION 



For 
has 



many years, 
been well 



the presence of intracluster (IC) stars 
documented by observa tions of dif- 
fuse light betw een cluster galaxies (e-K- |Zwicky | |l95l|: 
Gonzalez et alJl200ft IFeldmeier et, aDl2002 l : IFeldmeier et alJ 
20041: IZibetti fc White! 12004 iKrick et al. 112004) . Recently, 
numerous new searches for IC stars have uncovered individ- 
ual s ta rs between cluste r galaxi es (PNe: |Theuns fc Warrerl 
19971: lArnaboldi et alJ 120021: lArnaboldi et al. I 1200,4 
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IFeldmeier etafl Il99gl. 120031 : lOkamura et al. 1 120021: r ed gi- 
ants: iDurrell et alJl2002t SNela: lOal-Yam et al.11200.^. and 

indiv i dual tidal streams in c l usters jTr entham & Mobashor 
I l998t iGregg fc Wesd ll99St ICalcaneo-Roldan et alJ l200d : 
IFeldmeier et alJl2002l) . These studies, and others, all mea- 
sure an amount of intracluster light that is between 10% 
and 50% of the total cluster luminosity, suggesting that 
intracluster stars are a generic feature of galaxy clusters. 

The observable properties of intracluster stars are likely 
closely linked to the dynamical history of galaxy clusters, 
as IC stars are a natural by-product of cluster evolution 
in hierarchical models such as ACDM. Numerical simula- 
tions have demonstrated that the dynamical evolution of 
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galaxies in clusters is mainly dictated by global gravita- 
tional tides and fast encounters between galaxies that alter 



their morphologies and 


strip them of their dark matter ha- 


los (IMoore et al.l 


1996; 


lDubinskilll99Sl: iGhiena et alJll998l 


Calcaneo-Roldan et al. 
Die. IMoore et alJ (Il996l 


1200(4 iMaver et alJl200ll). For exam- 
) studied the dynamical evolution of 



individual galaxies inside a realistic cluster potential to show 
that, depending on their internal structure, they might un- 
dergo significant stripping of their stellar material. As galax- 
ies get stripped, stars in their tidal streams carry memory of 
the parent galaxy's orbit for several dynamical times. Spa- 
tial streams seen in Hydra and Co ma have thus been used 
to tra ce past dynamical interactions lICalcaneo-Roldan et alJ 
l200d) . 

As stripped stars become mixed into the global cluster 
potential, they still bear signatures of the cluster's accre- 
tion and dynamical history, although the extent of this is 
uncertain. Thus far, numerical studies that addressed such 
uncertainties have focused on the z = properties of in- 
tracluster stars. Dubinski (1998) used an approach where 
disk/halo galaxy models were substituted for dark matter 
halos in a cluster simulation in a critical universe to show 
that a significant fraction of stars end up unbound to in- 
dividual galaxies. Napolitano et al. (2003) compared a col- 
lisionless N-body simulation of the formation of a galaxy 
cluster with recent observations of Virgo IC PNe. A sim- 
ple scheme to identify a diffuse stellar component in the 
simulation showed that intracluster stars are largely unre- 
laxed in velocity and are clustered on 50 kpc scales at 500 
kpc from the cluster center, in good agreement with the 
observed clustering pr operties of the diffuse popula tion of 
PNe in Virgo. Recently. ISommer-Larsen et alJ (I2004S) looked 
at the z = photometric properties and metallicity distri- 
bution of IC stars in a cosmological simulat ion, and found 
them in good agreement with observations. IMurante et alJ 
( 2004) studied clusters in a large cosmological simulation 
and found fractions of u nbound cluster stars ra nging from 
10 - 50%. However, b oth ISommer-Larsen et al] <l2004f) and 
IMurante et al] i2004h 's spatial resolution was much coarser 
than the average scale size of cluster galaxies, which should 
overproduce intracluster stars, due to overly efficient strip- 
ping. 

The numerous remaining uncertainties in the expected 
properties of IC stars leave much room for improvement in 
this burgeoning field of study. For example, no consensus 
exists on when and how many IC stars should actually be 
produced in a case where the entire population of a forming 
galaxy cluster is realistically simulated. Nor is there consen- 
sus on whether IC stars originate mostly from the stripped 
outer parts of bright galaxies or rather from a population 
of tidally destroyed low surface brightness galaxies. It is 
also unclear whether observed cluster to cluster variations in 
measured IC fractions and radial profiles are due to differ- 
ent cluster formation histories and dynamical states, rather 
than a cosmic scatter in the expected properties of IC stars. 

Furthermore, little has been done theoretically to pre- 
dict the dynamics of the intracluster population; there is 
only one published analysis of the velocity distribution of IC 
PNe (ISommer-Larsen et al.ll2004h . Positions and velocities 
have been me asured for hundreds of IC planetary nebulae 
in Virgo (e.g. iFeldmeier et~al"ll2003l : lArnaboldi et aljl2002ft 
making it possible to use PNe to measure cluster mass, as 
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ICPNe outnumber cluster galaxies jFord et al.ll2002l) . How- 
ever, it is unclear if an accurate mass measurement with PNe 
is feasible, due to the substructure found in observations and 
simulations of IC stellar populations. 

In this paper, we address the outstanding questions 
above by using a simulation of a rich galaxy cluster in a 
realistic cosmological environment by using an unsurpassed 
coupling of high spatial and mass resolution. This allows us 
to follow the star formation histories of individual galaxies as 
faint as L*/10 for the first time, a crucial requirement to re- 
solve the formation of the large majority of the stars formed 
in the cluster. Nearly 90% of galaxy l ight in the Coma cluste r 
resides in galaxies with L> L./10 (Moba sher et al.ll2003h . 
Also, due to a spatial resolution higher than most previous 
analyses, our simulations produced a realistic distribution of 
galaxy sizes as a function of stellar mass and dark matter 
halo mass. This requirement is fundamental for a faithful 
description of the stripping processes in the dense cluster 
environment. Also, although simulations with a more ad 
hoc approach to tracing stars produce worthwhile results, 
a simulation that includes star formation produces a more 
accurate model of ICL as a continuous, ongoing process. 

The aims of this paper are: a) to shed light on the un- 
certain origin and properties of intracluster stars in a ACDM 
Universe and how they relate to the dynamical state of the 
cluster and b) test their use as dynamical tracers of a galaxy 
cluster potential. In §2, we describe the simulation; in §3, we 
describe the origin and evolution of the IC population; in §4, 
we discuss IC stars as tracers of the cluster's dynamical his- 
tory; and in §5, we evaluate the utility of IC PNe as tracers 
of cluster mass. 

2 THE SIMULATION 
2.1 GASOLINE 

We selected a candidate Coma cluster from an existing low 
resolution, dark matter only simulation in a concordance 
(A=0.7, f2o=0.3, (78=1) cosmology and resimulated it at 
higher resolution using th e volume renormalization tech- 
nique llKatz fc W hite 19 9j|). To perform this simulation, we 
used GASOLINE (TWadslev et al.ll2004l) . a smooth particle 
hydrodynamic (SPH), parallel treecode that is both spatially 
and temporally adaptive. The version of GASOLINE we 
used implemented i) Compton and radiative cooling, ii) star 
for mation and 'mini mal' supernova feedback as described 
by Katz et al. (1996), and iii) a UV background following 
lHaardt fc Madaul J1996J). The opening angle, 8, was 0.5 un- 
til z = 2 and 0.7 thereafter, and the tim e-stepping criterion, 
7/, was 0.2, as in lDiemand et ail i2004h . 

We ran four different simulations of the same cluster: 

- CI is a dark matter only run that we use to evaluate 
our halo completeness limit. This is a crucial measure as 
stars can form in poorly resolved halos but should not be 
included into our analysis. 

- C2 is our fiducial run that includes gas cooling pro- 
cesses, star formation and SN feedback. The high resolution 
region of this simulation contains dark and star particles of 
1.5 x 10 9 and 6.7 x 10 7 Mq respectively, with a force soft- 
ening of 3.75 kpc. There are 2 x 10 6 particles within the 
cluster's virial radius at z = 0. We will discuss results from 
C2, unless noted otherwise. 
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N dark 
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within R v i r 


Mq per particle 


kpc 


CI 


680814 




1.7 ■ 10 9 






3.75 


C2 


688754 


868995 


1.5 ■ 10 9 


7.2 


■ 10 7 


3.75 


C21ow 


85570 


143326 


1.2 ■ 10 10 


8.3 


• 10 8 


7.5 


C3 


679574 


516261 


1.5 ■ 10 9 


7.2 


■ 10 7 


3.75 



Table 1. Simulation Parameters 



- C21ow is a lower resolution version of C2 with 8 times 
fewer particles. 

- C3 is of equal resolution to C2, but includes the effects 
of additional feedback introduced by setting the gas at an 
entropy larger than 5 keV cm 2 at z = 3 to about 50 keV cm 2 , 
as measured by X-ray observations of clusters and groups 
(Borgani et al. 2002). 

All four clusters have R200, the radius enclosing an av- 
erage Sp/p = 200, of 2.2 Mpc, M200 ~ 1-2 ■ 10 15 M , and 
a ID velocity dispersion of 1000 km sec -1 . These are simi- 
lar to the observed properties of the Coma cluster jHughej 
Il989ft . 

Table |2~T1 gives the number of particles, particle masses, 
and force softenings of each simulation. The stellar masses 
given in the table are the maximum mass of a star within 
each z = cluster. 



2.2 Simulated Cluster Properties 

Figure shows the spatial distribution of the cluster stars 
at z = 0, represented by a random selection of 2% of all 
stars. The clear asymmetry in this distribution suggests that 
the galaxy population is unrelaxed, which we show in §3 to 
be the case, due to several major accretion events at low 
redshift. The outlined black squares show the size of a 0.5 x 
0.5 deg 2 field at 15 Mpc, the distance of Virgo. We discuss 
the properties of stars in these fiducial fields in §4. There are 
two large groups in the cluster: Group One is a Virgo sized 
group that entered the cluster 1-2 Gyr ago. Group Two is 
the size of Fornax, and entered the cluster 5 Gyr ago. 

20% of the simulated cluster's stars are unbound to in- 
dividual galaxies at z = 0, consistent with the values of 10 - 
50% suggested by observations (see §3.3 for additional dis- 
cussion). To distinguish stars gravitationally bound to indi- 
vidual galaxies from free-floating stars, we used the SKID 1 
halo finding algorithm on the gas and dark matter. We ini- 
tially ran SKID with a large linking length of 18 kpc, to 
minimize identifying loosely bound stars at z = as un- 
bound. Because that choice of linking length misses a few ha- 
los within the cD envelope, we also ran SKID with a smaller 
linking length of 9 kpc and included the additional galaxies 
in our sample. Based on a visual evaluation of the data, this 
combination produced the most robust identification of ha- 
los at z = 0. However, we found that our choice of linking 
length does not substantially affect the results. This proce- 
dure identified 650 galaxies within R,2oo- 
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Figure 1. The spatial distribution of stars in the simulated 
Coma-like galaxy cluster, C2. Only 2% of the star particles are 
displayed, for clarity. The outlined squares show the size of a 0.5 x 
0.5 deg 2 field for the cluster at a distance of 15 Mpc, the distance 
of Virgo. 



Due to the weak implementation of feedback in the 
simulations, star formation in dense gas regions proceeds 
unimpeded, resulting in a stellar baryon fraction of 36%, 
several times higher than t he 6 - 12% suggested by ob- 
servations of galaxy cl usters l|Bell et alJl2003t iBaloeh et alJ 
l200ltlCole et al. 1200 ill This "cooling crisis" is a well known 
weakness of the current generation of cosmological SPH 
simulations (Cole et al. 2001, Bell et al. 2003, Balogh 
et al. 2001, Kay 2001). Producing the correct amount of 
stars in a cosmological simulation is a problem that de- 
pends critically on numerical resolution and on the poorly 
known details of star formation and feedback from SNe 
(T hacker fc Couchman 2000, Bo rgani et al. 2001), although 
see ISpringel fc HernauisH i2003T> and ISommer-Larsen et ail 
( 2004}) for recent improvements. 

However, this overproduction of stars does not adversely 
affect our results because star particles in C2 galaxies have a 
realistic distribution, and are thus accurate tracers of the dy- 
namically stripped stellar population. If the simulated stel- 
lar halos are too compact or too diffuse, then the stripping 
efficiency by dynamical processes in the cluster would be 
under or overestimated, respectively. To determine the stel- 
lar distributions of galaxies in C2, we fit de Vaucouleurs 
effective radii (Re) to the 389 cluster galaxies with stellar 
mass > Miimit (10 9,9 Mq). Mu m it is the typical simulated 
stellar mass contained in a halo just above the conservative 
resolution limit of the simulation. We determine the halo 
resolution limit by a comparison with CI, the dark matter 
only run. The circular velocity function of CI started to flat- 
ten ~ v c = 100 km sec -1 , although C2's did not flatten until 
70 km sec -1 . Halos less massive than 100 km sec -1 would 
suffer less efficient cooling and may be incomplete due to 
numerical effects. We thus selected v c = 100 km sec -1 as 
a conservative resolution limit to be confident of complete- 
ness in our analysis. Nearly all galaxies with stellar masses 
greater than 10 9 ' 9 Mq are hosted by halos with v c > 100 km 
sec -1 . Halos more massive than this limit contain 98% of 
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Figure 2. Histogram comparison of the distribution of de Vau- 
couleurs effective radii for our simulated galaxies vs the distri- 
bution measured by Bernardi et al. (2003) for 9000 early type 
galaxies in the SDSS. 

the simulated stellar mass. We cut on stellar mass, rather 
than velocity or dark matter mass, to facilitate matching 
each simulated halo to a galaxy drawn from a Coma cluster 
luminosity function (§2.2.1). 

Figure |2] shows the resulting distribution of the sim- 
ulated galaxies' effe ctive radii, overp l otted with the distri- 
bution measured bv lBernardi et al~1 (J2003) , for 9000 early 
type galaxies in the SDSS. The R e distribution of C2 galax- 
ies is indeed in the same ran ge and peaks near the same 
scale size as that observed bv lBernardi et al~l (j^)03). The 
simulated galaxies are systematically slightly more compact 
than those observed, which may result in a lower unbound 
star fraction than reality (see §3.3). However, our measured 
distribution of scale sizes is a substantial improvement over 
existing studies of IC stars. We thus conclude that our simu- 
lated galaxies have sufficiently realistic stellar profiles to be 
used to study the production of intracluster stars in a rich 
cluster. 



2.2.1 The Galaxy Luminosity Function 

One of our main goals is to study what fraction of ICL orig- 
inates from galaxies of different masses. Our simulation re- 
solved a range of 4 decades in stellar mass; to determine 
how far down the galaxy luminosity function (LF) this stel- 
lar mass range includes and to evaluate which of our sim- 
ulations to focus this study on, we compare the luminosity 
functions of our simulations with that observed for the Coma 
cluster. Figure |3] shows a comparison between the original 
LFs of C2, C21ow, C3, and the lu minosity function for th e 
entire Coma cluster as observed bv lMobasher et all |2003). 
To convert the total stellar mass of each si mulated clus- 
ter ga laxy to an R magnitude, we combine the lSecker et alJ 
( 1997) obs erved R vs. B-R relat ionship for galaxies in Coma 
with the iBell fc de Jond fcOOlI) (M/L)_r vs ( B-R) relation 
with a scaled Salpeter IMF and ICole et ail j200Cfl galaxy 
evolution model. 

As star particles formed in the simulation act as dy- 
namical tracers, it is necessary that our simulation has a 
galaxy LF as close as possible to that observed for clus- 



Figure 3. Original luminosity functions of simulated cluster 
galaxies in three different runs: low resolution, high resolution, 
and high resolution with an entropy injection of 0.5 keV at z = 
3. For reference, the observed luminosity function of Mobasher et 
al (2003) is overplotted. 



ter galaxies. Figure [3] shows that C2 galaxies approximately 
follow a Schecter function, although the brightest galaxies 
overproduce stars. However, the lower resolution run, C21ow, 
neither produces the correct number of stars in massive ha- 
los nor resolves dark matter halos smaller than those that 
host galaxies. A run that resolves such a small fraction 
of the luminosity function would provide a very unrealistic 
and incomplete description of intracluster star accumula- 
tion. Although C3, the run with 0.5 keV particle -1 energy 
injection, produces a reasonable luminosity function above 
M|j, the formation of galaxies fainter than MJj is unrealisti- 
cally suppressed. This feed back recipe, w hile useful to study 
the properties of the IGM feorgani et al.l2 002). is too crude 
to accurately simulate galaxy formation within the cluster. 
Based on these comparisons, we choose to focus our analysis 
only on C2. 

To make our results directly comparable with obser- 
vations, we derive corrected stellar masses and R magni- 
tudes for each cluster g alaxy by normalizing to the observed 
iMobasher et al] (120031) luminosity function (LF) . Although 
the mass-luminosity relation of C2 is not correct, the basic 
correlation is correct: more stars form in more massive halos. 
Since the correlation between stellar and halo mass is mono- 
tonic, and we have a large number of halos, we are justified 
in using the observed luminosity function of the Coma clus- 
ter to calibrate the results of our star formation algorithm. 
We perform the normalization based on the 389 galaxies 
above the halo resolut ion limit (see above). The integrated 
IMobasher et al.l l|2003F) LF of the entire Coma cluster con- 
tains 389 galaxies brighter than R = -19.0, which includes 
86% of the total light they observed in cluster galaxies. We 
thus create rank ordered populations of 389 galaxies drawn 
from a Mobasher LF cutoff at R = -19.0, which is 2.8 mag- 
nitudes fainter than M^ of Coma galaxies. We convert their 
Mj?s t o estimated stellar masses by combining ISecker et al] 
( 1997) 's observed color-m agnitude relation for Coma with 
the lBell fc de Jond l)200ll^ color- (M/L) relation used above. 
We then fit a function relating these stellar masses to the 
rank ordered simulation stellar masses. We convert M_r to 
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Normalized Simulation LF 
Observed Mobasher et al LF 




Figure 4. The luminosity function of simulated cluster galax- 
ies after normalization, as compared to the observed luminosity 
function of Mobasher et al (2003). 



stellar mass, rather than just matching simulated stellar 
mass to M_r, so that we can normalize the mass of each indi- 
vidual star in the simulation. Figure0]shows the normalized 
luminosity function of our cluster. After normalizing, 10% of 
cluster baryon mass is in stars, vs. 36% pre-normalization. 

We used the stellar mass normalization derived at z = 
for cluster galaxies at all redshifts up to z = 1.1, correcting 
the derived luminosities for passive evolution. 



3 ORIGIN AND EVOLUTION OF THE 
INTRACLUSTER STARS 

Because this cosmological simulation includes a large dy- 
namic range of galaxies, we can derive the contribution of 
galaxies to the IC population, as a function of their abso- 
lute magnitude. In this section, we determine the origin of 
IC stars and use our result to predict and interpret observed 
colors and metallicities of IC stars. We also determine the 
evolution of the intracluster star fraction and evaluate the 
extent to which it is a function of the cluster's dynamical 
state. 



3.1 Determining the Progenitor Halos of the 
Intracluster Stars 

We traced all of the stars unbound at z = back until 
the epoch when they were last bound to a galaxy. As 85% 
of unbound stars were stripped since z = 1.1, we chose to 
only trace to that epoch. To reduce errors in assigning stars 
to their progenitor halos, we (i) do not allow stars to be 
"stripped" from the cD galaxy, but rather traced them back 
to their original halos and (m) require that a star is bound 
to its progenitor halo for at least two consecutive simulation 
outputs. We normalize the mass of each stripped star by 
the same factor that its progenitor's mass was normalized, 
as determined by the technique described in §2.2.1. We use 
these normalized masses in the subsequent results presented 
in this paper. 



3.2 Determining the z = Luminosities of the 
Progenitor Halos 

Thus far, observational studies of IC starlight have been lim- 
ited to relatively nearby, low redshift, clusters. Therefore, we 
determine the z — luminosity of the galaxy each progenitor 
evolves into. We trace the member stars of each progenitor 
galaxy forward in time to determine the halos they reside in 
at z = 0. We define the z = galaxy, G/, of each progenitor, 
Gi, as the z = galaxy that contains the highest fraction of 
Gi's stars. If the majority of Gi's stars are unbound at z = 
0, then Gf must contain > 40% of Gi's 10 most bound stars. 
If not, the progenitor halo is considered to be destroyed by 
z = 0. Only 1% of unbound stars from galaxies brighter 
than Mr — -19, were stripped from galaxies that were ulti- 
mately completely disrupted, and 9% of unbound stars were 
stripped from galaxies that ultimately merged with the cD. 

3.3 Evolution 

Observations of diffuse intracluster light have suggested 
that the majority of the IC star population in some clus- 
ters was accumulated early in the formation of the cluster 
jGonzalez et all200Cb . Authors have also suggested that the 
fraction of unbound stars in galax y clusters may be linke d to 
the dynamical state of the cluster jFeldmeier et al.l2002T) . In 
this section, we compare these past results to the evolution 
and fraction of intracluster stars in C2. 

Figure |S] shows the evolution of C2's intracluster popu- 
lation back to z = 1.1. The top panel compares the growth 
rate of the intracluster stars with the overall cluster growth 
rate. The peaks at z = 0.55 and z = 0.2 to 0.1 correspond 
to the infall of large groups, ranging in size from Fornax to 
Virgo. The growth of the IC population and the cluster as a 
whole are clearly linked and primarily occur by an ongoing 
series of events, rather than by a steady process. Stars are 
added to the intracluster halo both via stripping processes 
within the cluster environment and via the infall of large 
galaxy groups that already contain unbound stars. 

The bottom panel of Figure|K|shows the fraction of clus- 
ter stars that are unbound as a function of time. The IC star 
fraction in C2 displays a trend to increase with time. This 
trend is due to the continuous stripping of stars from clus- 
ter galaxies as they orbit. Although the infall of large galaxy 
groups adds stars to the IC population, the actual IC frac- 
tion is not necessarily linked to the dynamical state of the 
cluster. For example, the cluster underwent a substantial ac- 
cretion event at z = .55, but did not have an enhanced IC 
fraction as a result. This event increased the mass of the 
cluster's IC population primarily via stars already unbound 
within the infalling group, rather than by stripping within 
the cluster. The ratio of added IC star mass to added to- 
tal cluster mass was thus not larger than the ratio already 
present in the cluster, and thus did not increase the unbound 
star fraction. 

Since z = 1.1, the IC fraction varied between a mini- 
mum of 10% and a maximum of 22%, which closely matches 
the range in IC star fractions estimated for Virgo and sev- 
eral Abell clusters jFeldmeier et alJll998t iFeldmeier et alJ 
l2004h . although is less than the fraction estimated for 
Coma l| Bernstein et alJ Il995f) and a rich Abell cluster 
jFeldmeier et alJl2002jT This difference could be due to the 
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Figure 5. Top panel: Time evolution of the stripping and cluster 
growth rates compared with the time evolution of the fraction of 
cluster stars in the ICM. The values of stellar stripping and to- 
tal ICL growth are both normalized to the maximum ICL growth 
rate. The values of cluster growth rate are normalized to the max- 
imum cluster growth rate. Bottom panel: The time evolution of 
the unbound star fraction. 

facts that (1) the simulated galaxies are a bit more com- 
pact than those observed, and (2) the fraction we measure 
represents a lower limit to the fraction of cluster stars that 
are unbound, because can not include stars stripped from 
galaxies fainter than -19 in our analysis. 

Although stars are continuously added to the intraclus- 
ter population, the simulated cluster accumulated 50% of its 
IC stars within the last 2 Gyr. The upper panel of Figure 
03 shows that this recent increase in the number of IC stars 
was connected to a substantial accretion event experienced 
by the cluster between 1 and 2 Gyr ago. The mean dynam- 
ical time of this cluster is (3-7T /16Gp)~ 1//2 ~ 2 Gyr, thus the 
cluster is not in a dynamically relaxed state, as suggested 
by the asymmetric spatial distribution of the cluster seen in 
Figure 1. The variation of IC fraction with time in just one 
cluster suggests that there is no universal IC fraction. 



3.4 Origin 

Figure|S|shows the fraction of IC stars stripped as a function 
of the absolute magnitude of their progenitor halo, both at 
the time of stripping, z str ip, and at z = 0. 9% of stars are 
stripped from galaxies that become fainter than -19 by z 
= 0. The cumulative fraction of stars bound to halos with 
M_r < R at z = is also overplotted. The fact that the 
fraction of stars from halos with M_r < R at the time of 
stripping closely follows the fraction of stars contained in 
halos (plus an offset in Mr) shows that stripping efficiency is 
roughly constant with stellar mass. The offset in Mr results 
from our accounting for passive evolution when determining 
M_h at z > 0. The z = progenitor M_r distribution has 
a different shape than the z = z str ip distribution because 



stripped galaxies that are less massive than M* often merge 
into more massive galaxies. 

We find that massive galaxies contribute substantially 
to the IC population. Figure |S| also shows that 50% of free- 
floating cluster stars were stripped from progenitor galaxies 
with luminosities of M R or brighter. Because stars are pref- 
erentially stripped from the outer regions of galaxies, IC 
stars should bear the properties of stars in the outskirts of 
galaxies from which they were stripped. Assuming passive 
evolution, these should be similar to the z = properties of 
stars in the outskirts of progenitor galaxies. 

Stars in the outskirts of L* and brighter galaxies have 
photometric properties that differ from those of the stars 
producing the major ity of galaxy light, due to negative 
metallicity gradients l|Peletier et alJll99ff l. Therefore, stars 
stripped from the outer regions of the most massive galaxies 
could have the same colors and metallicities of an 'average' 
star in an intermediate luminosity galaxy. For example, the 
outer regions of elliptical galaxies have been observed to have 
(B-R) t hat are ^0.1-0.15 magnitudes bluer than the inner 
regions l|Peletier et alJll990h . lciadders et all ljl998h 's obser- 
vations of the red sequences of galaxy clusters showed that 
the average (B-R) of dwarf galaxies are ~ 0.3 magnitudes 
bluer than the average (B-R) of the most luminous cluster 
galaxies. Therefore, we expect many IC stars to have colors 
in common with intermediate luminosity galaxies (redder 
than the dwarfs, bluer than the most luminous), although 
they have been stripped from galaxie s brighter than M* . 
This prediction is consistent with the iDurrell et alJ l|2002f) 
measurement of -0.8 < [Fe/H] < -0.2 for the range of metal- 
licities of intracluster stars in Virgo. This metallicity range 
agrees with the mean metallicities of intermediate luminos- 
ity galaxies. 

Low surface brightness (LSB) and dwarf galaxies fall 
below the resolution limit of our simulation. We thus likely 
underestimate the fraction of unbound stars. Could LSBs 
and dwarfs ac tually be the popula tion contributing most of 
the IC stars? IDurrell et al. (2002) presented the example, 
that although dwarf galaxies contribute only ~ 10% of the 
light of Virgo, the IC stars contain ~ 20% of the light in 
Virgo. This means that for the IC stars to have all come 
from dwarfs, the number of dwarfs must have once been 3 
times their present number. This extreme scenario is not re- 
quired, as our results show that massive galaxies produce 
substantial ICL. Simulations that include galaxies with a 
more realistic distribution of profiles and/or include galax- 
ies even further down the luminosity function are necessary 
to know the stripping efficiency as a function of luminosity 
more precisely than we present here. 



4 INTRACLUSTER STARS AS TRACERS OF 
DYNAMICAL HISTORY 

In §3, we showed that the cluster is dynamically unrelaxed 
and has accumulated a substantial fraction of its IC popu- 
lation over the last few Gyrs. In this section, we investigate 
the extent to which this dynamical history is reflected in the 
observable spatial and velocity distributions of the cluster's 
stars. 
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Figure 6. The fraction of intracluster stars stripped from halos 
as a function of their progenitor's absolute magnitude, both at 
the time of stripping and traced forward to z = 0. The fraction 
of stars bound to halos with < R is also overplotted. 
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4.1 Spatial Distribution 

lOregg West! £M) and ICalcaneo-Roldan et al.l JUIJ) 
identified cold stellar streams in Coma with fj,R ~ 26.0. 
Such streams are evidence for recent stripping events, be- 
cause stars have not had sufficient time to get mixed into the 
global cluster potential. C2 underwent several major strip- 
ping events within the last 2 Gyr, associated with the infall 
of a Virgo sized group ('Group One' in Figure 1). As a re- 
sult, the cluster exhibits a wide stream of unbound stars 
that extends for hundreds of kpc between the cD galaxy 
and Group One, along the line that it recently crossed the 
cD. The average surface brightness of the unbound stars in 
this large feature is ~ — 26.5. We find no other streams 
in C2, but we likely underestimate the prevalence of streams 
that result from such an accretion event, because the galax- 
ies in this simulation are all spheroidal. Therefore, they all 
produce streams with internal velocity dispersions {(J stream) 
that are comparable to v ro t, rather than the much colder 
streams produced by disk galaxies (a stream ~ cr<Msfc). 

The surface brightness profile of C2's IC stars is in 
agreement with that observed in previous studies of both 
Coma and of other Abell clusters. Figure |7| displays the az- 
imuthally averaged surface density profiles of: the galaxy 
cluster (all stars+dark matter), the stars bound to individ- 
ual galaxies, and the intracluster stars. The profiles have 
been normalized to the average surface density of the clus- 
ter within the central 100 kpc. Although the recent merger 
event results in an asymmetric IC star distribution, their 
azimuthally averaged surface density profile beyond 250 kpc 
follows that of the gal axies, as also see n in tw o of the Abell 
clusters observed by iFeldmeier et alJ J2004f) . Within 250 
kpc, IC stars have a flatter profile than boun d stars (includ- 
ing cD stars) , at difference with the results of lMurante et all 
( 2004). In fact, our simulated IC population has a nearly flat 
surface brightness profiles between 50 and 250 kpc. Beyond 
50 kpc, most of the light in the bound stellar profile is not 
from cD stars. The average R-band surface brightness of 



Figure 7. Comparison of the surface density profiles of: the stars 
bound to galaxies, the intracluster stars, and all stars + dark mat- 
ter. These projected profiles only include stars within the virial 
radius of the galaxy cluster. 

intracluster stars 250 kpc from C2's cente r is 26.2, which 
agrees with that measured for Coma by iBernstein et alJ 
(1995). 

Figure|H|shows the azimuthally averaged total projected 
stellar profile of the cD galaxy with 2 de Vaucouleurs fits 
overplotted. The profile is inconsistent with a typical cD 
profile of a single R 1//4 law, either with or without an ex- 
cess at large radii. The unusually steep slope in the inner 25 
kpc is likely due to the unrealistic cooling flow and ongoing 
star formation at the center of our cD, and is not seen in the 
profiles of the other simulation galaxies. Although the radius 
ranges fit are arbitrary, the fits highlight a true flattening in 
the outer envelope of the cD galaxy, around 200 kpc. This 
flattening can not directly be due to Group One or Two, as 
they are a Mpc from the cluster cent er. This excess above a 
i? 1 / 4 profile is similar to that seen bv lFeldmeier et al l (120021) 
for the outer profile (> 80 kpc) of the cD galaxy in the rich 
galaxy cluster A1413. There is building evidence that such 
envelopes may even be ubiquitous features of cD galaxies of 
rich clusters (Gonzalez et al., in preparation). The fact that 
C2's cD display s an e xcess over R 1//4 at large R supports 
IFeldmeier et alJ ll2002h 's suggestion that such envelopes are 
a signature of recent cluster growth. An intracluster pop- 
ulation formed primarily during the initial cluster collapse 
differs from that of an intracluster population accumulated 
slowly over time via tidal processes, and is expected to fol- 
low an R 1 ' 4 out to even larg er radii than seen around C2 
jMerritdll98l |Pubinskilll998n . 

4.2 Substructure in Velocity Space 

We looked at the phase-space distribution of "planetary 
nebulae" (PNe) in four 120 x 120 kpc fields (30' x 30' at 
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Figure 8. Surface density profile of the cD galaxy. Two R 1 / 4 fits 
are overplotted; one was fit from 5-20 kpc and one was fit from 
40 - 200 kpc. 

a distance of 15 Mpc), selected by eye to lie outside the 
cluster galaxies. We used the luminosity-specific Planetary 
Nebula number dens ity of ai.o = 9.4 x 10~ 9 PNe Lg 1 
l|Ciardullo et, al]|l989T ) to randomly select a number of stars 
in each field representative of the number of Planetary Neb- 
ulae that existing spectroscopic surveys would observe. 

Figure|£j]shows the observable phase space distributions 
and line-of-sight velocity histograms of PNe in the two of the 
four fields plotted on Figure^ The fact that these 'observed' 
fields display both visible phase space correlations and non- 
Gaussian distributions of line of sight velocities reflects the 
young dynamical age o f the IC population, also found by 
iNapolitano et ahl (l2003f) when they applied the same anal- 
ysis to their simulations. The small scale phase-space sub- 
structure seen in our cluster's stars represents a lower limit 
on the actual substructure one would expect to observe, be- 
cause all of the galaxies in our cluster are elliptical, so stellar 
streams have higher internal velocity dispersions, on aver- 
age, and are thus more easily mixed than those that would 
result from a population that contains disk galaxies. 

4.3 The Anisotropy Profile of the IC Stars 

The top panel of Figure 1101 shows the ID radial profiles of 
the line-of-sight (los) velocity dispersions of IC stars, clus- 
ter galaxies, and cluster dark matter particles. To derive 
the error bars, we did a bootstrap resampling of the data 
for numerous cluster orientations. We find that the intra- 
cluster stars have a velocity dispersion on average ~ 20% 
smaller than that of either the galaxies or the dark mat- 
ter. The IC stars' velocity dispersion decreases with radius 
a bit more than that of the dark matter or the galaxies. 
This difference in velocity profile is due to the fact that, 
beyond 750 kpc, the intracluster star orbits are more ra- 
dial than those of either the dark matter or the galaxies. 
This bias is seen in the comparison between the anisotropy 
profiles (/? = 1— < erf > / < a, 2 >) of IC stars, galaxies, 
and dark matter in the lowe r panel of Figure 1101 Although 
ISommer-Larsen et al.l i2004h did not find such an anisotropy 
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Figure 10. Top panel: The line-of-sight velocity dispersion pro- 
files of the intracluster stars, the dark matter, and the cluster 
galaxies. /3 = — oo corresponds to circular orbits, = corresponds 
to isotropic orbits, and = 1 corresponds to radial orbits. Bottom 
panel: The 3D anisotropy profile. Cluster projection does not af- 
fect the 3D /3 profile. Error bars are not plotted for dark matter 
and unbound stars, because the formal bootstrapped errors are 
tiny, due to the large number of data points in each bin. 

bias in their simulations, it is not unexpected for unbound 
stars to have more radial orbits than galaxies or dark mat- 
ter; stars are preferentially stripped from galaxies on radial 
orbits because, on average, such orbits have sm aller percen- 
ters than isotropic orbits jTaffoni et alj feorm This bias in 
the anisotropy of intracluster stars, relative to that of either 
cluster galaxies or the dark matter, needs to be considered 
when constructing dynamical models based on intracluster 
stellar tracers. In the next section, we will do a simple cal- 
culation to demonstrate the possible effect of this difference 
on cluster mass estimates. 



5 INTRACLUSTER STARS AS A TRACER OF 
CLUSTER MASS 

The number of intracluster PNe easily surpass that of clus- 
ter galaxies, which suggests that PNe may be a more ef- 
fective tracer of t he cluster mass distribution than galaxies 
iFord et al .112002^) . Modeling the mass of a system based on 
the observed light and line of sight velocity dispersion pro- 
file relies on assumptions about either the mass distribution 
or the distribution of orbits. Therefore, the extent to which 
free floating stars can be used as mass tracers is determined 
by the extent to which we understand these properties of 
the intracluster population relative to that of the underly- 
ing mass distribution. 

In §4.2, we showed that IC stars are correlated in phase 
space and have a mean orbital distribution that differs from 
that of the cluster galaxies and dark matter. Furthermore, 
samples of galaxy cluster PNe are identified in discrete fields, 
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Figure 9. The spatial and velocity distribution of stars in the 2 'empty' fields in Figure^ located 500 kpc from the cluster center . 



such a s those of lArnaboldi et all i2002T) and lFeldmeier et alJ 
(1998). In this section, we use a crude approach to explore 
if PNe may nonetheless be useful tracers of cluster mass. 

The "projected virial theorem" is one simple dynami- 
cal method typicall y used to derive galaxy cluster masses 
iGirardi et alJ ll998l. However, numerous close pairs of ob- 
jects will be present in samples of Planetary Nebulae iden- 
tified in discrete fields, resulting in a projected virial mass 
that is more than an order of magnitude less than the ac- 
tual cluster mass. We therefore use the "projected mass" 
to evaluate Planetary Nebulae as ma ss tracers, because it is 
less sensitive to close pairs of objects llRines et alJl2003l) . We 
computed th e proje cted mass (PM) with the relation from 
iHeisler et al.l (|l985l) : 

i=l 

that assu mes isotropic orbi ts and a continuous mass dis- 
tribution iRines et alJl20o3l . Ri is each object's projected 
distance from the cluster center, Vi is each objects's observed 
line-of-sight velocity, and v is the average line-of-sight veloc- 
ity. The factor of 32 becomes 16 or 64 for circular or radial 
orbits, respectively. 

First, we calculated the projected mass with PNe from 
a number of 30' x 30' fields ranging from 1 - 14, and with 
a spectroscopic depth of either the brightest 0.5 or 1.0 mag- 
nitudes of the PNLF. We randomly placed 40 fields at dis- 
tances between 0.5 and 1.0 Mpc from the cluster center and 
automatically selected those fields farthest from any galaxy. 
The projected mass was thus calculated for numerous clus- 
ter projections. The average mass as a function of number 
of fields is shown in the left panel of Figure 1111 This plot 
shows that the number of PNe fields used in the analysis 
is more important than the depth of the spectroscopic fol- 
low up, for purposes of calculating mass. The scatter due to 
cluster projection is ~25% of the total cluster mass, similar 
to the dispersion found by iGovernato et alJ <l200ll) . Vary- 
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Figure 11. Left panel: Projected mass of the galaxy cluster cal- 
culated using a simulated sample of Planetary Nebulae selected 
in discrete fields between a projected radius of 500 and 1000 kpc, 
using either the brightest 0.5 mag or the brightest 1.0 mag in each 
of up to 14 fields. Right panel: Projected mass calculated using 
the N brightest galaxies, N randomly selected unbound stars, or 
N simulated Planetary Nebulae distributed between 14 0.5 x 0.5 
degree 2 fields. 

ing the clustercentric radii of the PNe fields showed that 
fields at a wide range of radii need to be included to avoid 
over/underestimates on the order of 20% the mean calcu- 
lated value. Coma is larger than Virgo , so it may be that 
PNe in fields su ch as lArnaboldi et all J2002lr s RCN1 and 
fields 2 and 6 of lFeldrrreler^t^lTT2OT3lj are at sufficiently 
large radii to be included in an accurate mass measurement. 

The projected mass calculated with 5 or more fields 
nearly asymptotes to the value calculated using unbound 
stars not restricted to individual fields (1.4- 10 1 Mq), shown 
in the right panel of Figure HT1 Surprisingly, this calculated 
PM is within 10% of the actual cluster mass of 1.5-10 15 . This 
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result is initially surprising because Figure 10 shows both 
that IC stars have a velocity dispersion that is somewhat 
smaller than that of the underlying mass distribution and 
that IC stars do not tend to follow isotropic orbits. These 
factors produce a mass underestimate relative to an isotropic 
population with the same velocity dispersion as the underly- 
ing mass distribution, when using Equation 1. However, this 
underestimate is canceled out by the fact that Equation 1 
systematically produces an overestimate of mass relative to 
true cluster mass, similar to the overestimate produced by 
the projected virial theorem without a surfa ce term correc- 
tion ijAceves fc Pereall999UR.ines et al.l2003h . We found this 
same overestimate when we calculated the average projected 
mass with C2 cluster galaxies, as seen in the right panel of 
Figure ITT1 The la scatter due to cluster projection is ~ 15% 
for the masses calculated using galaxies. We also found that 
the projected masses calculated using galaxies are very sim- 
ilar to the projected virial masses without a correction for 
the surface term. 

Our analysis shows that IC PNe are worth a closer in- 
vestigation into their utility in making detailed mass models 
of nearby clusters. The fact that our analysis yielded a clus- 
ter mass very close to the true mass was partially a fortuitous 
canceling of two effects; a detailed analysis would thus likely 
produce a more reliable approach to cluster mass modeling 
with individual IC stars than we have presented here. 



6 CONCLUSION 

In this paper, we used a high-resolution, SPH simulation of a 
Coma-like cluster formed in a cosmological context to study 
the formation, evolution and properties of intracluster stars 
in a rich galaxy cluster. This simulation resolves galaxies 
as faint as Mr = —19, which includes over 85% of the light 
from cluster galaxies. Furthermore, the simulated stellar dis- 
tributions have a range of effective radii that matches that 
observed by Bernardi et al. (2003). This dynamic range and 
accurate stellar distribution are essential for relying on sim- 
ulated stars as accurate tracers of the dynamically stripped 
stellar population. 

Overall, we find that the fraction and distribution of IC 
stars in a hierarchical scenario are a good match to the ob- 
served properties, qualitatively simil ar to the results found 
in t he lower resolution simulation s of lMurante et all <l2004l) 
and ISommer-Larsen et all J2004J) . ~ 20% of stars are un- 
bound to any galaxy at z = 0. The ICL accumulation is an 
ubiquitous, ongoing process, with an unbound stars fraction 
that slowly increases with time. The intracluster population 
is formed via both stripping within the cluster and via the 
infall of large groups that contain stars that are already un- 
bound. As a result, there is a link between the accumulation 
of IC stars and infall into the cluster. However, the actual 
IC star fraction at any given time is not necessarily corre- 
lated with the dynamical state of the cluster, with both the 
lowest and highest IC fractions occurring after substantial 
merging events. 

Our analysis traced the fraction of stars stripped from 
galaxies as a function of galaxy luminosity. We found that 
the most luminous galaxies in a cluster contribute a sub- 
stantial fraction of stars to the unbound population, and 
that the stripping efficiency is roughly constant over a wide 



range of galaxy masses. As a result, we predict that IC stars 
will have photometric properties in common with the av- 
erage properties of intermediate luminosity galaxies, which 
is consistent with existing observations (e.g. Durrell et al 
2002). 

We verified the phase space and velocity distribu- 
tion of IC stars found in previous numerical work. IC 
stars exhibit non-Gaussian velocity distributions, as well 
as visible correlations in their observable phase space plots 
jNapolitano et al]l2003ft . We also showed that the orbits of 
IC stars are more radially biased than those of the overall 
galaxy or dark matter, which results in the IC stars having a 
lower l.o.s. velocity dispersion in the outskirts of the cluster 
than in its inner regions. 

We find that despite the fact that IC stars do ex- 
hibit phase space substructure, and have l.o.s. velocities and 
anisotropies that differ from that of the underlying distri- 
bution, they are promising tracers of galaxy cluster mass. 
The depth of the spectroscopic observations in each field 
makes little difference in the calculated masses. The pro- 
jected masses calculated with discrete fields of unbound stars 
quickly approaches that calculated with randomly selected 
unbound stars, which is within 10% of the true cluster mass 
due to a fortuitous canceling of two effects. Our analysis 
suggests that more detailed analysis may provide a "correc- 
tion" to the projected mass found with IC stars, similar to 
the known correction for the projected virial mass estimate. 

Much progress remains to be done to complete our un- 
derstanding of intracluster stars and what they can tell us 
about the formation history of the clusters which they in- 
habit. Some groups are already (e.g. Mihos et al 2004) start- 
ing to include galaxies with more realistic profiles into nu- 
merical simulations, which will enable a more quantitative 
analysis of substructure predictions and of stripping efficien- 
cies. Better resolution coupled with a more realistic imple- 
mentation of feedback should be the goal of future cosmolog- 
ical simulations used to study the properties of intracluster 
stars. These improvements will enable more robust predic- 
tions of colors, gradients and light profiles than possible with 
the current generation of simulations. 
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